Coulomb drag in graphene: perturbation theory 
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We study the effect of Coulomb drag between two closely positioned graphene monolayers. In 
the limit of weak electron-electron interaction and small inter-layer spacing (/ii(2),T ^ '"/d) the 
drag is described by a universal function of the chemical potentials of the layers Mi(2) measured in 
the units of temperature T. When both layers are tuned close to the Dirac point, then the drag 
coefficient is proportional to the product of the chemical potentials po oc l-i.ifJ'2- In the opposite limit 
of low temperature the drag is inversely proportional to both chemical potentials pn oc T^ /{pifi2)- 
In the mixed case where the chemical potentials of the two layers belong to the opposite limits 
A^i ^ T <C /i2 we find po oc pi/p.2- For stronger interaction and larger values of d the drag 
coefficient acquires logarithmic corrections and can no longer be described by a power law. Further 
logarithmic corrections are due to the energy dependence of the impurity scattering time in graphene 
(for pi{2) ^ T these are small and may be neglected). In the case of strongly doped (or gated) 
graphene pi{2) ^ v / d ^ T the drag coefficient acquires additional dependence on the inter-layer 
spacing and we recover the usual Fermi-liquid result if the screening length is smaller than d. 

PACS numbers: 72.80.Vp, 73.23.Ad, 73.63.Bd 
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Transport measurements are conceptually simplest and 
by far the most common experimental tools for studying 
inner workings of solids. Within linear response, the out- 
come of such measurements is determined by the proper- 
ties of the unperturbed system, which are often the ob- 
ject of study. In a typical experiment a current is driven 
through a conductor and the voltage drop along the con- 
ductor is measured. In most conventional conductors at 
low temperatures the resulting resistance is mostly deter- 
mined by disorder (which is always present in any sam- 
ple), while interactions between charge carriers lead to 
corrections that affect the temperature dependence^. 

Consider now a drag measurement in a bi-layer system 
consisting of two closely spaced but electrically isolated 
conductors^. Passing a current Ji through one of these 
conductors ( "the active layer" ) is known to induce a volt- 
age V2 in the other conductor ("the passive layer"). The 
ratio of this voltage to the driving current pjj = — V2//1 
(the transresistivity or the drag coefficient) is a measure 
of inter-layer interaction. At low enough temperatures 
the drag is dominated by the direct Coulomb interaction 
between charge carriers in both layers^. 

The physics of the Coulomb drag is well understood if 
both layers are in the Fermi liquid state'^''*. The current 
in the passive layer is created by exciting electron-hole 
pairs (each pair consisting of an occupied state above the 
Fermi surface and an empty state below) in a state char- 
acterized by finite momentum. The momentum comes 
from the electron-hole excitations in the active layer cre- 
ated by the driving current. The momentum transfer is 
due to the inter-layer Coulomb interaction. Therefore it 
follows from the usual phase-space considerations that 
the drag coefficient is proportional to the square of the 



temperature oc T^. Remarkably, this simple argu- 
ment is sufficient to describe the observed temperature 
dependence of pu (deviations from the quadratic depen- 
dence are due to the effect of phonons)^. 

The phase-space argument however does not describe 
the physics of the effect completely. Indeed, in the pas- 
sive layer the momentum is transferred equally to elec- 
trons and holes so that the resulting state can carry cur- 
rent only in the case of electron-hole asymmetry. Like- 
wise, this asymmetry is necessary for the current-carrying 
state in the active layer to be characterized by non- 
zero total momentum. In the Fermi-liquid theory the 
electron-hole asymmetry can be expressed^ (assuming 
either a constant impurity scattering time or diffusive 
transport^) as a derivative of the single-layer conductiv- 
ity (Ji[2) with respect to the chemical potential. In con- 
ventional semiconductors^, the asymmetry appears due 
to curvature of the conduction band spectrum (leading 
to the energy dependence of the density of states (DoS) 
and/or diffusion coefficient): dai/djii sa aij pi. 

Theoretical calculations^*^ typically focus on the drag 
conductivity ctu. The experimentally measurable drag 
coefficient pn is then obtained by inverting the 2x2 
conductivity matrix. To the lowest order in the inter- 
layer interaction'^ (assuming cr/) ^ <^\(2)) obtains 



PD = o'_D/(a-icr2). 



(1) 



Combining the above arguments (and assuming that the 
single-layer conductivities Oi are given by the Drudc for- 
mula) we arrive at the Fermi-liquid result^ 



PD 



e2 ^1^2 



-A, 



(2) 
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where A12 depends on the matrix elements of the inter- 
layer interaction, the Fermi momenta of the two layers, 
and the inter-layer spacing d (in the diffusive regime, 
where the mean-free path £ <^ d, A12 contains an ad- 
ditional logarithmic dependences.). The precise form of 
A12 is well known and can be obtained by means of either 
the diagrammatic formalism^ or the kinetic equationi^>^. 

Recently drag measurements were performed in a sys- 
tem of two parallel graphene sheets^. It was shown that 
this system offers much greater flexibility compared to 
prior experiments in semiconductor heterostructuresi^. 
The drag coefhcient depends on the following parameters: 
(i) temperature T, (ii) chemical potentials of the lay- 
ers /Xi, (iii) inter-layer spacing d, (iv) mean- free path 4, 
and (v) the interaction strength. Earlier experiments.'-'^ 
were performed on samples with large inter-layer spacing, 
such that /i 3> w/d (where v is the Fermi velocity). In 
contrast, the graphene-based system allows one to scan 
a wide range of chemical potentials (by electrostatically 
controlling carrier density) from the Fermi-liquid regime 
with v/d,T to the Dirac point = where the 

drag vanishes due to electron-hole symmetry. 

Most experiments^iiS (including that of Ref.'o') are per- 
formed in the ballistic regime £i 3> rf, where po does not 
explicitly depend on disorder, albeit the conductivities 
cr_D and ai do. The graphene-based sample of Ref. H is 
characterized by much smaller inter-layer spacing in com- 
parison to previous experiments^ , with all data taken at 

T < v/d. (3) 

Therefore in this paper we do not consider temperatures 
larger than the inverse inter-layer spacing, even though 
our approach remains valid for T > v/d. 

In this paper we address the problem of the Coulomb 
drag in graphene by means of the perturbation theory. 
In the leading order, i.e. in the limit of weak interaction 

a = e^/v < 1, (4) 

the drag conductivity ao is described by the standard 
(Aslamazov-Larkin-type, see Fig. [T]) diagram^. Unlike 
usual metals, the single-layer conductivity in graphene 
comprises two competing contributions: one due to dis- 
order and one due to electron-electron interactioii^i^— . 
We assume that the dominant scattering mechanism is 
due to disorder. Hence throughout the paper we assume 

a^Tr < K Tr, (5) 

where r is the impurity scattering time (in the case of 
energy-dependent impurity scattering time this parame- 
ter should be understood as r(niax[/i, T]), see Sec lIVBl) . 
The latter inequality ensures that the system is in the 
ballistic regime. Under this condition^i the single-layer 
conductivities cr 1(^2) are again given by the Drude-like for- 
mula with the impurity scattering time. As we will show, 
<^d/<^i{2) ^ a^Tr <^ 1, which allows us to evaluate the 
drag coefhcient pn using Eq. ([T]). Our theory can be 
further extended to the case of stronger interaction or 



TABLE I: Asymptotic expressions for the drag coefBcient (|6]) 
in the limit of weak inter-layer interaction and small inter- 
layer spacing ^1^2), T <C v/d. 

parameter region drag coefBcient 

^il,P2<.T PC ^ 1.41 a^-^^TT^ Eq.JMI 

T«Mi<M2 PD^a^^^^lnt^ Eq.m 



vanishing disorder, where the condition ([5]) can be lifted. 
However it is more convenient to perform such calcula- 
tions within the framework of the kinetic equatio n^^'^^ . 
The results of that work are reported elsewherei^. 

In the limit of small inter-layer spacing pi(2)^ T v/d 
the resulting drag coefficient is the "universal" function 

PD = o?^ ro{pi/T,p2/T). (6) 

The function rQ[xi,X2) is characterized by the limiting 
cases Pi ^ T and pi ^ T, see Table HI In the vicinity 
of the Dirac point pu oc pip2/T^, while in the opposite 
limit p£) cx T"^ / {P1P2) (with additional logarithmic fac- 
tors, see Sec. IIII|) . In the "mixed" case /zi <C T <C ^2 
we find p£) cx pi/p2- This regime is apparently realized 
in the experiment of Ref. when the bottom layer is in 
proximity of the Dirac point. 

For relatively large inter-layer spacing, pi ^ v/d (and 
thus Pi ^ T), a new regime appears, where the drag coef- 
ficient acquires additional dependence on the inter-layer 
spacing. Here the assumption ^ may be relaxed as the 
perturbation theory remains valid also for intermediate 
values of a (see Sec. IIV Al where we assume pi = p2). 
StiU, for ^ < r we find po oc p^ /T'^, a,t p ^ T the drag 
coefficient reaches its maximum, and further decays for 
p^ T. The latter regime is described by a long crossover 
from the above logarithmic behavior to the Fermi-liquid 
result which is only achieved in the limit of the small 
screening length xd 3> 1. Thus pnip ^ ^1'^) cannot be 
described by a single power law. 

The remainder of this paper is organized as follows. In 
Sec. U we present the perturbative calculation that gives 
the drag resistivity in the limit a — ^ and — > 0. While 
the resulting expression can only be evaluated numeri- 
cally, we analyze all interesting limits analytically for the 
case of two identical layers in Sec. [ill and also for the ex- 
perimentally relevant case where the two layers are char- 
acterized by different chemical potentials, see Sec. IIIII 
In Sec. IIVI we discuss drag at intermediate values of a 
and d as well as the role of energy-dependent scattering 
time. The paper is concluded by a brief summary and 
a discussion of the experimental relevance of our results. 
Throughout the paper we use the units with h = 1 and 
only restore the Planck's constant in the results for pc- 
Technical details are relegated to Appendices. 
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passive layer (2) 



active layer (1) 



FIG. 1: [Color online] The lowest-order diagram for drag 
conductivity. Solid lines represent the exact single-electrons 
Green's functions in the presence of disorder. Spatial coordi- 
nates are labeled by the numbers shown at the corners of the 
triangles and correspond to the subscripts in Eq. (|13|l . Wavy 
lines represent the inter-layer interaction. 



B. Non-linear susceptibility in graphene 

The non-linear susceptibility of electrons r^ (a;) is a 
response function relating a voltage y(ri)e*"* to a dc cur- 
rent it induces by the quadratic response: 

1 = Jdr.j dr2T{uj;r,,r^)Vir,)V{r2), (11) 

with I being the induced dc current. From gauge invari- 
ance it follows that 

j driT{uj]ri,r2) = J dr2T{uj;ri,r2) = 0. 



PERTURBATIVE CALCULATION OF THE 
COULOMB DRAG 



In terms of exact Green's functions of a disordered con- 
ductor the non-linear susceptibility can be written as^^i^ 



Consider the limit of weak interaction a — ?■ 0. In this 
case the drag conductivity ao can be calculated with the 
help of the lowest-order diagram? 



a/3 



IGttT ^ 
q 



du! 



sinh 2y 



R |2 



(7) 

where the subscripts 1 and 2 refer to the active and pas- 
sive layers respectively, is the retarded propagator of 
the inter-layer interaction, and r"(w, q) is the non-linear 
susceptibility (or the rectification coefficient). 



A. Inter-layer interaction for small a and d 

The bare Coulomb potential has the usual form 
27re2 _ 2Tre^ 



Vii = V2: 



Vv. 



-qd 



q 



(8) 



The dynamically screened (within RPA) inter-layer prop- 
agator can be written asi^ 



1 



nf 



(9) 

where 11^ is the single-layer retarded polarization opera- 
tor. For Of ^ 1 screening is ineffective and the inter-layer 
interaction is essentially unscreened. Moreover, in the 
limit d — we may disregard the exponential and write 
the interaction propagator as 



^12 — 



(10) 



As we will show below, the non-linear susceptibility T de- 
cays exponentially for q 3> max(r, /i). Therefore, in view 
of Eq. ^ the limit d — > is equivalent to the condition 
ix<^v/d (the case > v/d is discussed in Sec. IIV Ap . 



27r 



tanh 



tanh 



2T 



tanh | 'yi2\^\^) 



2T 



tanh 



e — cj — /i 
2T 



2T 



72i(e;-w) 



,(12) 



where the arguments of T are suppressed for brevity, the 
subscripts 1 and 2 refer to the spatial coordinates ri and 
r2 , respectively, and the triangular vertex 7 is given by 



7i2(e;w)= Gf2(e + cj)-G^^2(e + ^) 



(J9 



^23 (e)J3G^i(e). 

(13) 

Here J3 is the current operator and integration over the 
coordinate is assumed (see Fig. [1]). 

The traditional approach calls for averaging the vertex 
7 over disorder. Within the Fermi-liquid theory the av- 
eraged 7 does not depend on ei^. This result is based on 
the usual approximation that the most important contri- 
bution comes from electrons near the Fermi surface. In 
contrast, in graphene- based systems the regime where the 
chemical potential is smaller that temperature fii <^ T 
is accessible. Then the vertex 7 retains its dependence 
on e and should be evaluated with carei^r— . Details of 
the calculation are provided in Appendices [K\ and [B] The 
result is 



■j{e,uj,q) = —2qevT{e 



2„2 



uj'^ + 2euj — v^q 



(14) 



/[2e + a;]'-«V 

^2g2 _ ^2 sgn(e + w)0o(e,^,g), 



where the fimction O^^e, uj, q) ensures that the expression 
under the square root in Eq. (1141) is always positive (the 
minus signs in 0o(e,a;,(7) reflect the contribution of both 
electrons and holes and appear after summing over all 
branches of the Dirac spectrum in graphene) 
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9o(e,a;,(?) = B{vq - \2e + uj\) [e{-LO - vq) - 9{uj - vq)] + 9{vq - \uj\) [9{2e + uj-vq) - 9{-vq - 2e - uj)] . (15) 



Approxiniating^i the impurity scattering time by a constant r (the effect of energy-dependent T(e) is discussed below 
in Sec. IIVB[) . we can express the non-hnear susceptibility in terms of dimensionless variables 



and find the following form 



W Q 



>,q) = -2^q.9(M^,Q;|), 



(16a) 



^ I J^,,S.^,, l2{z;W,Q;x), \W\ > Q 



Z\/i2~~T 



I{z;W,Q;x), \W\ < Q 



(16b) 



/(z; W, Q; x) — tanh Q ^ tanh ^ 1_ tanh — ^ ^ — — tanh — ^ — ^ (16c) 



From Eq. (|16cp it is clear that at the Dirac point the 
non-linear susceptibility vanishes: 



I{x^O) = => Till = 0) = 0. 



(17) 



Thus there is no drag at the Dirac point (physically, due 
to electron- hole symmetry). 



C. Perturbative results for the Coulomb drag in 
graphene 

1. Drag conductivity in graphene 

Using Eqs. (|10l) and (|16l) we can find the drag conduc- 
tivity in graphene ([7]). As usual for an isotropic system 
in the absence of external magnetic fields = 5°'^cjd- 
Now, in the limit a — > and d — > we find 

ao^aeTTjay—.—j, (18a) 

where the dimensionless function fo{xi,X2) is defined as 



/o(a;i,a;2) = ~2 j J 



dW 



sinh W 



9ix,)gix2). (18b) 



where we have suppressed the arguments of the function 
g {W,Q;x) for brevity. 

In general the function fo{xi, X2) has to be computed 
numerically. Below, we evaluate this function analyti- 
cally in the physically interesting limiting cases of large 
X (the low temperature regime) and small x (the vicinity 
of the Dirac point). 



2. Smgle-layer conductivity in graphene 

In order to find the drag coefficient po one needs to 
know the single-layer conductivity ai . Under our assump- 
tions the single-layer conductivity is completely deter- 
mined by the weak impurity scattering and can be writ- 
ten in the form 



where 



00 
2 f 
ho(x) = — dz 
71- J 



^Trh, 



\z\ 



o(f) 



(19) 



cosh 



2^ I X, X > 1, 
TT I 21n2, a; < 1. 



(20) 



3. The drag coefficient 



Finally, using Eqs. (UJ, (|18ap . and (HH) we find that 
the drag coefficient can indeed be expressed in the form 
with the dimensionless function ro(a;i,X2) defined as 



ro{xi,X2) 



foixi,X2) 
ho{xi)ho{x2) 



(21) 



In the case of two identical layers this function depends 
on one variable only 



''0(2:) 

and is shown in Fig. [2l 



fo{x,x) 
hlix) ■ 
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FIG. 2: [Color online] The dimensionless drag coefficient 
ro{x,x) (two identical layers). The solid line shows the re- 
sult of numerical evaluation using Eqs. (|21|l . (I20p . (|18b[l . and 
(|16|l . The red dashed line shows the asymptotic behavior (|27|) 
in the vicinity of the Dirac point. The green short-dashed line 
shows the asymptotic behavior (|33p at small temperatures. 



The expressions ^ and ([2T|) give the perturbative re- 
sult for the drag coefficient in graphene-based bi-layer 
systems in terms of the dimensionless functions ho{x) 
and fo{xi,X2). The applicability of this "universal" re- 
sult and its experimental relevance is discussed in Sec. lIVI 



II. 



ASYMPTOTIC BEHAVIOR OF THE DRAG 
BETWEEN TWO IDENTICAL LAYERS 



Consider now the case where the two layers in the drag 
experiment are identical, i.e. are kept at the same tem- 
perature and chemical potential. This case has not been 
yet realized experimentally in graphene-based systems, 
but has a long history of theoretical researcb^i^. 



A. Vicinity of the Dirac point 

We start with situation where the graphene sheets are 
tuned close to the Dirac point, i.e. 

Given that at the Dirac point the nonlinear susceptibility 
vanishes [see Eq. ([T7|) ]. we can expand it to the lowest 
order in the small parameter /i/T. Expanding first the 
expression (I16cp we find 



I{x -> 0) « -4a;- 



sinh W sinh zQ 



(22) 



(cosh zQ + cosh W) 
This corresponds to the quadratic expansion of Eq. (jI8bp 

foix,x)^Mix^. (23) 



The numerical coefficient A^i can now be determined by 
using the approximation ([^ in Eqs. (|16bp and (|16b|) and 
then evaluating the integral in Eq. (jI8bp . The result is 



Ml « 1.1. 



(24) 



Therefore the drag conductivity in the vicinity of the 
Dirac point is independent of temperature and is propor- 
tional to the square of the chemical potential 



(25) 



In this case the drag conductivity is independent of T. 

The single-layer conductivity can be found by expand- 
ing the function ho{x). Thus in the vicinity of the Dirac 
point the conductivity of a single graphene sheet due to 
impurity scattering [see Eq. ([5])] iai^ 



41n2 



Therefore the drag coefficient is 



Pd(^ < T) « 



e2 161n2 2T2 



(26) 



(27) 



This result is represented in Fig.[2]by the dashed red line. 



B. Low temperature limit 



In the opposite limit 



we notice that the function I2 is given by Eq. (|22p with 
interchanged W and x and can be written as 



I{x > 1) 



m d 



sinh a; 



Q dz cosh zQ + cosh : 



(28) 



Then the integral over z in Eq. (jl6b|) is dominated by 
z 1 so that the algebraic function of z in the integrand 
may be approximated by unity. The corresponding con- 
tribution (jl6bp to the non-linear susceptibility may now 
be approximated by the expressioni^^ 



gix,\W\<Q)=A 



W 



•1- 



sinhx 



cosh Q + cosh x 



(29) 



Now the momentum integral in Eq. (|18bp is logarithmic 
and is dominated by large values of momentum Q ^ W. 
Then the function (jl8bp takes the form 



64 

fo{x,x) = 



dWW^ fdQ 



sinh^ x 



sinh W J Q (cosh Q + cosh : 



\2 ■ 



W 



(30) 

The ratio of the hyperbolic functions in Eq. (pO|) is similar 
to the step function: it's equal to unity for Q <^ x and 
vanishes at larger values of momentum Q ^ x. Therefore 
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X effectively acts as tlie upper cut-off and tfie momentum 
integral can be approximated by a logarithm 



dQ 



sinh 



w 



Q (cosh Q + cosh a;)^ 



In- 



ly 



Therefore, in the low temperature limit the leading con- 
tribution to the drag conductivity is logarithmic in the 
chemical potential and quadratic in temperature 



a,,(T«M)«f «VTVln^. 



(31) 



The single-layer conductivity at low temperatures is 
determined by the chemical potential^^ and thus is given 
by the Drude formula 



(To ~ (2/7r) e^/^T. 



(32) 



Consequently the drag coefficient is similar to Eq. ([2]). 



Y 7T m — . 

e2 3 T 



(33) 



This is to be expected, since at low temperatures T <C 
the phase-space argument yielding the dependence is 
justified and the electron- hole asymmetry determines the 
dependence on the chemical potential. The logarithmic 
factor is of course beyond such qualitative estimates. 

This result ([55)1 is represented in Fig. [5] by the dashed 
green line and is shown together with the result of the di- 
rect numerical calculation (shown by the solid blue line) . 
Note that the drag conductivity (I5T|) was calculated with 
logarithmic accuracy. 



III. ASYMPTOTIC BEHAVIOR OF THE DRAG 
BETWEEN INEQUIVALENT LAYERS 



The dimensionless function /i {x) is characterized by the 
two limits. If the second layer is also close to the Dirac 
point, then 

fi{x < 1) wMx, 

which yields a straightforward generalization of Eq. (|23p : 
/o(a;i,X2) ~ A/1X1X2. If, on the other hand the chem- 
ical potential of the second layer is large compared to 
T, then the drag conductivity (jl8b[) can be calculated 
by using the linear approximation (j22p for g{xi) and the 
low temperature approximation (|29p for g{x2). Due to 
the exponential decay of Eq. (f22|) at Q ^ W , the mo- 
mentum integral is dominated by Q ^ 1 and thus the 
combination of the hyperbolic functions in Eq. (|29p can 
be replaced by unity. Thus in this limit the drag con- 
ductivity in independent of the chemical potential of the 
second layer and 



/i(a;> 1) «AA2, AA2«3.26. 



(34) 



Now, the single-layer conductivity in the first layer is 
given by Eq. (|26p. while in the second layer we should 
consider both limits. If ^2 is also small then the resulting 
drag coefficient is a trivial generalization of Eq. ([27]) : 



/0I3(M1,M2 < T) 



e2 T2 



(35) 



In the opposite limit the drag conductivity is independent 
of the properties of the second layer, since the integrals 
in Eq. (jl8bp are dominated by the region where both fre- 
quency and momentum are of order T. The single-layer 
conductivity in the second layer however is determined 
by ^2 [see Eq. (|32|) ] and thus 

8 In 2 /i2 /i2 



/Od(M1 < T < /i2) 



Consider now the more realistic^ situation where the 
two layers are characterized by different chemical poten- 
tials. We will still assume that the temperatures of the 
two layers are the same. 

Note that in this Section the subscripts of the chemi- 
cal potentials /ii and /12 do not indicate the passive and 
active layers, but rather simply distinguish between two 
layers with different carrier density. 



A. One layer near the Dirac point 

Firstly, suppose that one of the two layers is charac- 
terized by a small chemical potential 



Ml < T. 



Then in Eq. (jl8bl) the function g{xi) may be expanded 
using Eq. 



B. One layer with high carrier density 

Secondly, if the chemical potential in one of the layers is 
much larger than temperature (without loss of generality 
we may also assume that the chemical potential of the 
other layer is also smaller) 

M2 > 7", M2 > Ml: 

then arguing along the lines of the previous subsection 
we find that the drag conductivity is independent of the 
largest chemical potential 

ao = a2e2TV/2(Mi/T). 

Again, we characterize the dimensionless function f2{x) 
by the two limits. The situation when the second layer is 
near its Dirac point was already discussed in the previous 
subsection. Therefore, similar to Eq. ([M)) 

f2{x < 1)kN2X. 
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The drag coefRcient in this case is given by Eq. ((36|) . 

In the case where > Mi ^ ^ the calculation is 
similar to that presented in Sec. Ill Bl However now the 
combination of the hyperbolic functions [coming from the 
approximation (j29p ] in the momentum integral similar to 
Eq. pop comprises two step functions and the integration 
is cut off by the smaller chemical potential. Hence 

/2(a:»l)«yln^, 

the drag conductivity is independent of the larger chemi- 
cal potential /Lt2 . The drag coefficient in the limit of large 
(but smaller than /ii) is a generalization of Eq. 



PD(/i2 > Ml > T) w a' 



— ^ In—. (37) 



rs 
rs 

53 




IV. BEYOND THE LOWEST ORDER 
PERTURBATION THEORY 

A. Finite inter-layer spacing and static screening 

Let us now discuss what happens for non-zero inter- 
layer spacing d. For simplicity, we will consider the case 
of two identical layers. Generalization to the case of in- 
equivalent layers is achieved along the lines of Sec. IHIl 



FIG. 3: [Color online] The sketch of the drag conductivity 
(in the units of a^e^r'^) as a function of the chemical po- 
tential illustrating the results of Section IIV A II in the case 
T <^ Nav/d. The blue line corresponds to the quadratic de- 
pendence (|25|l in the vicinity of the Dirac point. The four 
regions at large chemical potentials correspond to the four re- 
gions discussed in Section [IV A II If T 3> Nav/d, then the 
region (iia) should be replaced by (ii) , the logarithmic depen- 
dence ln(l/a) should be replaced by ln[v/{Td)], and the 
limits v/d and T/{Na) should be exchanged. 



1. Vanishing interaction strength 

In the limit a — >■ the inter-layer spacing appears only 
in the exponential factor in the unscreened inter-layer 
interaction ([5]). Therefore, the momentum integration 
acquires a firm upper cut-off q < v/d, but the behavior 
at small momenta is unchanged. 

As we have seen in Sec. HIl the behavior of the drag 
coefficient in the vicinity of the Dirac point is determined 
by frequencies and momenta of the order of T (within our 
perturbation theory the single-layer conductivities in the 
two layers are independent of each other and d). Under 
the assumption ^ these momenta are small compared 
to v/d: taking into account the exponential in Eq. ^ 
yields a small correction to the numerical coefficient in 
Eqs. (HZl) and (03: 1.41 ^ 1.41 + iMTd/v. 

In the opposite limit of low temperature we have found 
the behavior of the drag conductivity that is determined 
by large momenta. In the result (j331) the upper limit 
of the logarithm is given by the chemical potential [due 
to the step-like behavior of the non-linear susceptibility 
(j29p ]. while the lower limit is given by temperature as 
the typical value of frequency in Eq. ([50]). 

Clearly, if finite d is taken into account then the loga- 
rithmic behavior in Eq. p3p may change: the upper limit 
of the logarithm will now depend on the relative value of 
the chemical potential and inverse inter-layer spacing. 

At the same time the lower limit of the logarithm also 
depends on our approximations. Indeed, so far we have 



considered the effect of the unscreened inter-layer inter- 
action. Consider now the role of the static screening. If 
both d and a are small then we may approximate the 
inter-layer interaction ^ by the expression 



^12 



27re2 



2n'^ 



(38) 



At low temperatures the leading contribution to the 
static polarization operator is 



n> = o) = ^. 

TTV 



(39) 



Thus the interaction propagator can be written as 



^12 



2'Kav 



vq + 2Nafi 



-qd 



Here = 4 is due to spin and valley degeneracy. The 
unscreened interaction is a good approximation as long as 
the inverse screening length is small compared to typical 
momenta, i.e. as long as Aa/i ^ T. In the opposite 
regime the lower limit of the logarithm in Eq. (|3ip will 
be determined by the screening length. 

Combining the above arguments we conclude that in- 
creasing the chemical potential the following four regimes 
may be gradually achieved. 

(i) Na^ ^ T <C /U <C v/d. This regime was considered 
in Sec. [H] leading to Eqs. ^ and 
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(ii) Na^ <^ T <^ v/d <C /z. If the chemical poten- 
tial is increased beyond the inverse inter-layer spacing, 
then the momentum integration in Eq. (pO|) is cut off by 
v/d instead of /i. The logarithmic behavior of the drag 
conductivity Eq. ([31]) will be modified and cr^i no longer 
depends on the chemical potential 



a'e'T'r'ln—. 



(40) 



(iia) T <^ Nafj, ^ /i ^ v/d. Depending on the actual 
values of T, d, and a it is possible that Na/i exceeds tem- 
perature while the chemical potential is still smaller than 
the inverse inter-layer spacing /i ^ v/d. Then instead of 
the previous regime we find 



a^e^T^T^ In 



1 

Na' 



(41) 



(iii) T <C Nafi ^ v/d ^ fi. Increasing the chemi- 
cal potential further leads to the regime where the static 
screening can no longer be neglected. Now the lower inte- 
gration limit in Eq. (j30p is effectively given by the inverse 
screening length rather than the frequency. The upper 
limit is still determined by the inter-layer spacing. There- 
fore the drag conductivity again depends logarithmically 
on the chemical potential^^ 



CTD 



(42) 



but now this is a decreasing function, indicating the exis- 
tence of the absolute maximum of the drag conductivity 
as a function of the chemical potential. 

(iv) T < w/d < Nafi < Finally, if the chem- 
ical potential is so large that the screening length be- 
comes smaller than the inter-layer spacing the momen- 
tum integral in Eq. pop is no longer logarithmic. In 
this regime wc recover the standard Fermi-liquid result^. 
Note, that in this regime the step-like combination of the 
hyperbolic function in the non-linear susceptibility (j29p 
is completely ineffective and may be replaced by unity. 
Given that in this regime the integration is dominated by 
momenta large compared to temperature the non-linear 
susceptibility may be further linearized in frequency. The 
resulting expression 



C(3) eV^r^ 

4 (fcj.d)2(>fd)2' 



0-Z5 = 



>^ = Aakp, 



(43) 



differs from that of Ref . 3 only by the factor reflecting val- 
ley degeneracy in graphene. Such expression for the drag 
conductivity was previously obtained in Ref. [20l . The 
results of this subsection are illustrated in Fig |31 



2. Intermediate interaction strength and T 

The above results rely on the smallness of the inter- 
action strength a. However, if Na > 1, then (i) the 
approximation psp might not be justified and we would 



need to use the full expression ([9]) for the interaction 
propagator; (ii) the four regimes specified in the previ- 
ous subsection may not exist, since it might happen that 
T/{Na) < T < v/{Nad) < v/d. In this case there are 
only two distinctive regimes (for fi ^ T): (a) /i ^ w/d, 
and (b) fi'^ v/d. The latter regime is usually identified 
with the Fermi-hquid lesuli^^^^^^^i^^, i.e. Eq. 

In this subsection we derive the approximate expres- 
sion for the drag conductivity for large values of the 
chemical potential ^ » T, which accounts for the pos- 
sibility of Na > 1 (possible for small a <C 1, but large 
iV ^ 1; here we still consider identical layers). In this 
regime the single-layer conductivity (|32]) is large (since 
fiT >• Tt 1) and therefore we can still limit our con- 
sideration to the diagram in Fig.[TJ Moreover, for /i ^ T 
we can somewhat relax the condition ([5]) on the interac- 
tion strength and require the electron-electron scattering 
time to be larger than the impurity scattering time 



> T 



T > — 



2ti „ M 



Now we can follow the usual steps leading to the Fermi- 
liquid result (l43l) . We consider only the static screening 
by approximating the polarization operator by Eq. p9p 
(see Appendix [C] for details). We further assume^^ 
that the dominant contribution to the drag conductivity 
comes from the region vq > u. Finally, we assume that 
in that region the result of the momentum integral is de- 
termined by the upper integration limit and is therefore 
independent of uj. This allows us to evaluate the fre- 
quency integral and represent the drag conductivity in 
terms of the single integral over momenta. Under these 
assumptions we find similarly to Eqs. ([50)1 



Z 'Im'Z Z p It'*' Td \ 

(ID = a e T T /o — ;a; — , 



T' 



V J 



(44a) 



Jo{x;a;X) « — 



dQQ-^e 



3„-4AQ 



3 7 [{Q + a{x)Y - a{xYe-^^QX 
sinh^ X 



(cosh Q -|- cosh x) 



where 



a{x) 



N 



ax. 



(44b) 



Here we have approximated the function g{x, \W\ < Q) 
by Eq. ((29l) , but neglected the frequency under the square 
root. In the limit fj, ^ T the combination of the hy- 
perbolic functions in g{x, \W\ < Q) has the form of the 
step function, which effectively cuts off the integration at 
Q ^ X. Thus we arrive at the approximate expression 



fo{x;a;\) ~Y J 



dQQ^e 



3„-4AQ 



[{Q + a{x)Y - a{xYe-^^Q] 



(44c) 



9 



The results of the previous subsection can now be re- 
covered for a ^ 1 by neglecting the terms proportional 
to o? in the denominator. In contrast, the Fermi-liquid 
result (|33|) can be obtained by (i) assuming a ^ 1 and 
keeping only the terms proportional to o? in the denom- 
inator, (ii) assuming x ^ 1/4A and thus replacing the 
upper integration limit by infinity, and (iii) replacing the 
lower integration limit by zero. Moreover, for /i 3> w/c? 
the integration limit can be extended such that the result 
becomes a function of one single parameter 



/o(x;a; A) w /o(4Aci!), 



where 



/o(y) = -17 



3 7 [{Z + vf^f~e-^X 



(44d) 



(44e) 



The function (|44ep describes the crossover between the 
regimes (iii) and (iv) of the previous subsection (see 
Fig.[3|). This can be seen by evaluating the integral in the 
two limits (here 70 « 0.577216 is the Euler's constant) 



/o(j;«l) 



32 

y 



In- 



1 



7o 
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(45a) 



a=0.01, Td/v.O.l 



direct numerical evalu: 
static screening 
Eg. (44e) 




X = jU/ T 



FIG. 4: [Color online] Results of the numerical evaluation of 
the drag coefficient for a = 0.01 and Td/v = 0.1. The squares 
represent the calculation of Eq. ((Tj) with the only approxima- 
tion that the polarization operator in the screened inter-layer 
interaction ((9)1 was evaluated in the absence of disorder. The 
red line corresponds to the same calculation, but replacing 
the polarization operator by Eq. (|39p . i.e. taking into ac- 
count only static screening. The blue line corresponds to the 
approximate expression (|44e|l .valid for /j, ^ v/d. 



/o(.»l)-^. (45b) 

y 

It turns out that by numerical reasons this crossover 
spans a large interval of values of the chemical poten- 
tial such that the Fermi-liquid result (l43t is practically 
unattainable in graphene-based drag measurements®. We 
illustrate this point in Figs. H]- [S] 

First of all we evaluate the whole expression ([7]). In 
order to do that we need to evaluate the full non-linear 
susceptibility (jl6p . In addition we need to determine the 
polarization operator in graphene at finite temperature 
and chemical potential. In the ballistic regime we ne- 
glect the effect of disorder on the polarization operator. 
This is the only approximation in this calculation, see 
Appendix \C\ for details. The role of the disorder is dis- 
cussed in the following subsections. The corresponding 
data is shown in Figs. H]- [5] by blue squares. Further re- 
sults of the numerical evaluation of the drag conductivity 
([7]) are shown in Appendix |D] 

Then we can simplify calculations by using the static 
approximation p9|) for the polarization operator. The re- 
sult of this calculation is shown in Figs. |l]-[n]by the red 
line. For weak interaction a = 0.01 and relatively small 
inter-layer spacing (such that Td/v = 0.1) the results 
of the static screening approximation are indistinguish- 
able from the "full" calculation as can be seen in Fig. 2] 
At the same time, for "intermediate" values^ii^^ a = 0.3, 
Td/v = 0.2, the static screening approximation "works" 
well for small and large chemical potentials, while some- 
what overestimating the overall peak height. 



a = 0.3, Td Iv = 0.2 




1 1 — I 1 1 — I ' ' — 

1 10 100 1000 

x = /j/T 



FIG. 5: [Color online] Results of the numerical evaluation of 
the drag conductivity for a — 0.3 and Td/v = 0.2 shown 
in the log-log scale. The straight green line represents the 
Fermi- liquid result (|43l) . The solid blue line corresponds to 
the approximate expression (|44e[) . valid for jj, ^ v/d. 



Finally, we can evaluate the approximate expression 
(|44e[) . This expression was derived assuming /i ^ v/d. 
For the parameter values used in Figs. S] - [S] the ex- 
pression (|44el) fits the exact calculation for n/T > 10. 
However, the Fermi-liquid asymptotic (represented by 
the straight green line in Fig. [SJ is not reached until 
/i/T > 200. Clearly this is very far from the param- 
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eter regime relevant to the experimental observation of 
the Coulomb drag in graphene^;^. It is therefore not sur- 
prising that the Fermi-liquid-like approximations overes- 
timate the observed values of the drag^. The data in 
Fig. [5] show that the asymptotic results, such as Eq. p3)) 
may sometimes be achieved only at the extreme values 
of parameters. In order to describe the effect in the in- 
termediate (or realistic) parameter regime, one needs to 
evaluate Eq. ([7]) with only minimal approximations. 



In the simple case T(e) 
simplifies to 



const, the function K(z) 



K{z) 



-2z 



1 - T^VQ2 



and we recover Eqs. (fT6b| and (|T6b)) . 

2. Vicinity of the Dirac point 



(47) 



B. Energy-dependent scattering time 

The drag coefficient ([2T|l and the asymptotic results 
of Sec. ini and Sec. IIIH as well as the numerical data 
of Sec. IIV Al were obtained assuming a constant im- 
purity scattering time r. In graphene, the scattering 
time strongly depends on the type of the impurities and 
usually depends on energy^^"— . In the context of the 
Coulomb drag in graphene a similar issue (namely, the 
momentum-dependent scattering time) was investigated 
in Ref. [2^ in the framework of the kinetic equation. 



1. Non-linear susceptibility 

Consider now the effect of the energy (or momen- 
tum) dependence of the scattering time. Within our as- 
sumption ([5]) we can still consider the ballistic Green's 
functions (jAlOp . only now instead of being the overall 
factor, the scattering time is a part of the integrand 
in the non- linear susceptibility p^ . Due to the 6- 
function form of the Green's functions we can focus on 
the energy-dependent r regardless of the microscopic im- 
purity model. 

Repeating the steps leading to Eq. (IT6| with the 
energy-dependent r(e) we find that the functions (|16b|) 
and (|16b[) become (see Appendices \X\ and [Bl for details) 



/ 1 



dz- 



9{x) = < 



I{z)K{z), \W\>Q, 



(46a) 



Jz"^ - 1 

dz ^. -. I{z)K{z), \W\<Q, 



2./1- 



where 



K{z) 



zW -Qt{T[zQ-W]) 
zQ-W ^^(T) 



(46b) 



zW + Qt{--T[zQ + W]) 



zQ + W 



-{T) 



so that the explicit factor of the scattering time in 
Eq. (|16a|) should be understood as t(T). This choice 
of the prefactor is not essential and is dictated by the 
discussion of the case <C T below. 



Taking into account energy dependence of the scatter- 
ing time T(e) modifies the non-linear susceptibility and 
thus changes the Coulomb drag. In the limit fi -^T the 
region Q ^ W becomes important. As can be seen from 
Eqs. (|16bp and (|16b[) . if T(e) = const, then precisely at 
Q = |iy| the non-linear susceptibility vanishes. Other- 
wise, Eqs. (j46[) may contain a divergence. 

Consider for example Coulomb scatterers. Then2£iS 



Tie) 



(48) 



In this case the function (|T7)l depends on the relation 
between \W\ and Q 



K{z) 



Q, Q>W, 
z\Wl Q<W, 



(49) 



and as a result the integral (jl8bl) contains a logarithmic 
divergence at Q = since now 



(50) 



Same result holds for the case of strong short-ranged im- 
purities that also yield the linear dependence r ^ |e| (up 
to logarithmic renormalization which is inessential in the 
present context)^. 

Similarly, in the case of weak short-ranged disorder^-S 

r-i(e)=7k|, (51) 

the function (|46b[) does not vanish^^ for Q — \W\: 

z^Q^ + W^- 2z^W^ 



K{z) 



2 < 



Q- 



z\W\ 



[z^Q^ 
z^Q^ + 



- 2Q2 



(z2 



1^2 



>\Wl 



< \w\ 



Therefore the non-linear susceptibility contains the same 
logarithmic divergence (150]) . 

Finally, even in the case of logarithmic renormalization 
of the scattering time^i^ the function ([Tf]) does not van- 
ish for Q — \W\ leading to Eq. ([5D|) . We conclude, that 
Eq. (jSOl) is the generic behavior, while the vanishing (for 
Q — \ W\) non-linear susceptibility ([T6| is the artifact of 
the approximation r(e) = const. 

The divergence indicates that in the region W ^ Q 
the 5-function approximation for the Green's functions 
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(jAlOp is invalid. Going beyond this approximation is 
qualitatively equivalent to regularizing the divergence by 
the scattering time t{T). Thus the drag coefficient in 
the vicinity of the Dirac point acquires an additional log- 
arithmic factor 



2 ft A^lM2 
J^2 



lnTT{T). (52) 



3. Low temperature limit 

In the low temperature limit the non-linear susceptibil- 
ity is determined by electrons near the Fermi surface and 
one can set e in the triangular vertex (1131) to be equal 
to the chemical potential. Then it appears that in all 
subsequent expressions one has to replace r(e) by t(/x) 



r(e) w t(^). 



(53) 



Let us understand this statement in more detail. Con- 
sider for example Coulomb scatterers. In the limit T 
the non- linear susceptibility is dominated hy Q ^ W . In 
this case we can evaluate the integral in Eq. (|46a|) using 
the function K{z) from Eq. (j49p and approximating I{z) 
by Eq. (gHl)- As a result we find 



giQ » W) 



AW sinh : 



x+Q 



cosh Q + cosh X 



~Q l + e^+Q 



Oix-Q). 



(54) 



Therefore 



The momentum integral is still logarithmic, 
we can approximate g{Q ^ W) by 

AW 

giQ » W) « —X, 



which indeed recovers the statement ((5^ [notice that 
t{T)x = t{h)]. 

This argument is not restricted to Coulomb scatterers 
since in the limit x ^ Q the function I{z) in Eq. ([25)1 is 
essentially a (5-function: I{z) w —AW6{x — zQ). There- 
fore the function K{z) in the integrand in Eq. be- 
comes K{z) — > —2t{ii)/[zt{T)] for an arbitrary T(e). 

In addition, the function g contains the divergence (|50|) 
aX Q ^ W which is always present for T(e) ^ const. Sim- 
ilarly to Eq. (|5^ this results in an additional logarithmic 
contribution (|50p . which is typically subleading. For ex- 
ample, the Fermi-liciuid result (|43p acquires an additional 
factor 1 -I- {Td/v)\nTT{T). Thus we conclude that for 
/i 3> r taking into account the energy dependence of the 
impurity scattering time does not affect our results. 

Finally, the divergence ([50)) appears only if both chem- 
ical potentials are small. Indeed, if /ui <C T ^ ^2, then 
while the non-linear susceptibility in the first layer is 
given by Eq. (I46p . in the second layer the scattering time 
may be replaced by its value at the chemical potential 



and thus Eq. (jl6p applies. Now the divergent denomina- 
tor \W'^ — Q^l^^/^ of the former susceptibility is canceled 
by precisely the same factor in the numerator in the lat- 
ter. Hence, the product g{xi)g{x2) is finite at \W\ = Q. 



C. Plasmon contribution 

The above results were obtained while neglecting the 
plasmon pole of the interaction propagator. Let's now 
estimate the plasmon contribution keeping a small. Con- 
sider for simplicity identical layers. In the limit d — > 0, 
the plasmon pole in the propagator of the inter-layer in- 
teraction is very similar to that of the single layer. In- 
deed, introducing dimensionless functions 



- 

■^12 — 



27re2 



-D, n" = 



2q 



q 



(55) 



the inter-layer propagator ([9]) can be expressed as 
D-i ^ (1 + /3P)2e«'' - ;32p2g-gd^ /3=— . 

TT 

In the limit d — >■ the pole corresponds to the solution 
of the equation 

1 + 2/3P = 0, 

which differs from its single-layer counterpart by the fac- 
tor of 2 only. 

Now the drag conductivity (O is determined by the 
product of the square of the interaction propagator and 
two non-linear susceptibilities. Using the dimensionless 
notations (|55|) and (|16|) we can write the integrand in 
Eq. ([7]) in the form 



P\TlR |2 ^ „2_2„2 '/ H „2|£)| 



(X e T a 
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(56) 



Calculations of Sec. II CI are essentially equivalent to ar- 
guments based on the Fermi Golden Rule yielding the 
perturbative result po ~ 01^ ■ Now we show that taking 
into account the plasmon contribution results in addi- 
tional smallness, justifying the perturbative calculation 
of Sec. EC] 

Plasmon modes in a single graphene sheet in the vicin- 
ity of the Dirac point were studied in Ref. [ij. The plas- 
mon pole appears in the region 

W>Q. 

Adjusting for the above factor of 2, the plasmon disper- 
sion and decay rate are given by 



W„ 



Q 



Q + 26l 



-(g-l-a), a = 8aln2, (57) 



r _ Q 

^ 161n2 VQ + 2a 



3/2 



(Q + 5). 



(58) 
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The plasmon contribution to the drag conductivity 
may be described by the following form of the inter-layer 
interaction propagator 



\D\ 



1 



At the same time, as argued in Ref. Ilj], the typical mo- 
menta dominating the relaxation rates are not too small 



> 



a. 



(59) 



Then the decay rate is small Tp ^ and the interaction 
propagator has the form of a sharp peak. Estimating the 
typical frequency at the peak by the plasmon dispersion 
we note that 



Q + 2a 



(60) 



and therefore the inter-layer interaction takes the form 
(omitting inessential numerical factors) 



\D\ 



(W^ - Wp)2 -I- a 



Consider now the non-linear susceptibility g. Typi- 
cal momenta ([55]) are small enough, which allows us to 
approximate the non-linear susceptibility (jl6bp by its 
asymptotic value in the small Q — > limit 



5i(Q^0) 
As a result, for Q ~ a, 

Ad? (X 



(W^ - Wpf -f <56 



a-^S{W-Wp). 



(61) 



For larger momenta Q ~ 1 the result is similar. The 
interaction propagator is now given by 



\D[ 



{W - WpY + 



At the same time the non-linear susceptibility (jl6p con- 
tains a factor ^/W^ — Q"^ . Using Eq. (pOj) . we find 



[W ~ Wpf + (54 



<5M(W^- Wp). 



If however, the temperature dependence of the scat- 
tering time is taken into account, then instead of van- 
ishing at = Q the non- linear susceptibility contains 
the divergence (1501) and thus a more careful analysis is 
necessary. Now we need to find the plasmon dispersion 
and decay rate at T ^ in the presence of disorder. This 
complicated problem lies beyond the scope of the present 
paper. Here we estimate the plasmon contribution to ou 
in the most relevant region Q \W\ T. It turns out 



that under our assumption ([5]) this contribution is small 
compared to the leading approximation ([T5)) . 

In order to find the plasmon dispersion we need to cal- 
culate the polarization operator at T 7^ in the presence 
of disorder. This can be done with the help of the ki- 
netic equation derived in Ref. Isil . In comparison to the 
usual r-approximation this equation contains an addi- 
tional term describing the suppression of backscattering 
in graphene. It turns out that in the region Q ^ \W\ ~ T 
the polarization operator may be approximated by 



dn 



(62) 



where the thermodynamic density of states is given by 



dn 



- In 



2 cosh — 
2T 



(63) 



The absence of the extra term l/r in the denomina- 
tor (c.f. Ref. IH) is precisely due to the suppression of 
backscattering: for arbitrary q and uj there is indeed a 
rather involved expression generalizing this term to the 
case of graphene, but ior Q ^ \W\ ^ T this contribution 
is small and may be neglected. 

Solving for the plasmon dispersion with the help of 
Eq. ^ we find 



Wp^Q, 



rp = T- 



Thus the plasmon dispersion is changed little from 
Eq. ([57)) . where for Q ^ a we also have Wp « Q. How- 
ever the decay rate is now completely determined by dis- 
order. Now the interaction propagator takes the form 



\D\ 



l/Ti 



Q2 {w -WpY + i/t-^t'^ 



Multiplying this expression by the diverging , we notice 
that corrections to the linear plasmon dispersion are still 
determined by the interaction as in Eq. (15(11) . Therefore, 
we find that the resulting contribution contains the small 
factor o^Tt and is negligible under our assumption ([5]). 

In order to estimate the drag conductivity we now need 
to integrate the product g2|£)|2 over frequency and mo- 
mentum, see Eqs. ((T]) and (j56l) . This product contains 
a small factor of at least ce^, or, if energy dependence 
of the scattering time is taken into account, a small pa- 
rameter o^Tt. Therefore the plasmons contribute in the 
subleading order in the perturbative expansion in a. 

The conclusions of this Section are confirmed by com- 
paring the results of numerical evaluation oi po using ei- 
ther the full dynamically screened interaction ([9]) or only 
the static screening, see Eq. ([M]). The results are illus- 
trated in Figs, mo the difference between the two results 
is only noticable for larger values of a. Thus taking into 
account plasmons does not lead to any new qualitative 
features of the theory. 
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D. Spectrum renormalization 

If Coulomb interaction is taken into aecount, then 
the Dirac spectrum in graphene acquires logarithmic 
corrections^^. This can be understood in terms of the 
renormalization of the interaction parameter and 
disorder strcngth'^^. The renormalization group flow ter- 
minates at max(/z, T) and at lower energy scales we can 
treat the parameters a and v as scale-independent and 
equal to their renormalized values. The disorder scat- 
tering time retains its explicit energy dependence which 
follows from the microscopic impurity model. 

In our calculation of the drag conductivity all fre- 
quency integrals are effectively cut off by temperature, 
while the momentum integrals are cut off by either T or 
fj,, whichever is larger. In all of these cases we can treat 
the spectrum as linear with the renormalized velocity. 
Then our result (|52|) is still applicable, with the veloc- 
ity and interaction parameter a taking the renormalized 
values f[max(/i, r)] and a[max(/i, T)]. 



E. Experimental relevance 

1. Carrier density 

Experimental results^ are expressed as a function of 
carrier density rather than the chemical potential as we 
have discussed in this paper. The relation between the 
carrier density n and the chemical potential fi can be 
obtained by integrating the density of states p(e): 



J dep{e) [riF(e; ^Ji) ~ ^^(e; 0)] , 



(64) 



where np{e\ /i) is the Fermi distribution function. 
In graphene p{e) = 2|e|/(7rw^) and the integral 



de\e\ 



e e — /i 

tanh tanh 

2T 2T 



(65) 



can be easily evaluated in the limiting cases [cf. Eq. (1551) ] 



Tiv^ \ (41n2)^r, M < 



(66) 



which of course recovers the T = expression for /i ^ T. 

If impurity scattering is taken into account then the 
density of states in the vicinity of the Dirac point satu- 
rates to a value determined by disorder— 



n(/.,T<r-i) = -^. 



(67) 



However for Tt 3> 1 this effect is not important. 

In experiment the carrier density may be obtained from 
measurements of the Hall coefficient in a non-quantizing 



magnetic field H. In graphene this is more compli- 
cated than in usual metals since the Hall coefficient van- 
ishes at the Dirac poinli^ due to electron-hole symmetry. 
While at low temperatures the behavior of the Hall co- 
efficient is rather complicated'^^, at high temperatures 
we can use the conventional Boltzmann kinetic equa- 
tion with the energy- dependent cyclotron frequency^ 
a;c(e) = eHv'^ /(ce). Then we find 



2 2 

-e V 



2 2 

~e V 



de 



de 



dnpie) p{e)Ttr{e) 
de l + ojj{e)Tl{ey 



(68a) 



(68b) 



Using Eqs. (|68p at low temperatures fj,'3> T and for weak 
magnetic fields H we of course recover the classic result 



1 



^ly)H nec 



Rh 



with the electron density given by Eq. (|66|) . 

Exactly at the Dirac point the Hall conductivity van- 
ishes, (Jxy[p- = 0) = 0, as can be seen directly from 
Eq. (I68bp : all functions in the integrand, except for a;c(e) 
are even in e. For finite p, the Hall coefficient is lin- 
ear in the chemical potential and thus linear in carrier 
density, as can be seen from Eq. (j66p . 



Rh fx 



pv 
ecT^ 



The numerical coefficient in the above expression de- 
pends in the precise nature of impurities (see Sec. IIVB|) . 



2. Single-gate setup 

At the time of writing, there is only one published re- 
port of a Coulomb drag measurement in graphene-based 
double-layer system^. In this experiment there is only 
one gate controlling the carrier density in both layers. 
The carrier densities can then be found by solving two 
electro-static equations^i^ 

eVsG = Ms + e^(«s + "-t)/Ci, (69a) 

PB= PT + e^nT/C2, (69b) 

where Vbg is the voltage applied to the bottom gate, ps 
and pt are the chemical potentials of the bottom and 
top layer respectively, Ci is the capacitance of the oxide 
layer between the gate and the bottom layer, and C2 is 
the capacitance of the inter-layer spacing. 

Eqs. were used in Ref. |25( to deduce that the gate 
voltage Vbg is proportional to the carrier density. Let us 
estimate the density for which electrical and chemical 
potentials of the bottom layer become comparable: 



(70) 
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Here the numerical value is estimated using the param- 
eters of the experimental devic o^'^^ . Thus the linear re- 
lation Vbg oc is valid for ub ^ n*p , which is satisfied 
for all densities considered in Ref. [2^ 

In the vicinity of the Dirac point and in the presence 
of disorder [see Eq. (j67|) ] all terms in Eqs. (|69|) are lin- 
ear in carrier density. As a result, both ub and nx are 
proportional to the gate voltage Vbg and in particular 

Vbg oc ^iB■ 

According to Ref. the carrier density in the top layer 
depends on Vbg only weakly and remains finite when the 
bottom layer is tuned to the vicinity of the Dirac point, 
such that fiB < ^J'T■ In such conditions the drag coeffi- 
cient is described by Eq. (1551) . Since in this regime the 
gate voltage seems to affect mostly the bottom gate of the 
device used in Ref. @, we conclude that when the bottom 
layer is tuned towards the Dirac point, the drag should 
vanish linearly with the gate voltage. This conclusion is 
consistent with the experimental results of Ref. 

At the large carrier densities the drag should vanish 
as some power of the gate voltage as shown in Fig. [2] 
Since in the experiment the Fermi wavelength and the 
inter-layer spacing are of the same order of magnitude, 
the decay of the drag coefficient at large fi is described 
by Eq. ([37|. Neglecting the weak dependence of ut on 
the gate voltage reported in Ref. we conclude that 

PD ~ 

which qualitatively agrees with the experimental results. 

Given that the drag vanishes both at small and large 
fi (or gate voltages, or carrier densities), there must be 
a maximum at some intermediate value of fj,, which is 
determined by temperature and sample geometry. Thus 
the results of the perturbation theory, as shown in Fig. [2l 
qualitatively describe all features of the drag observed in 
the experiment. Moreover, for relatively small values of 
the interaction parameter^^ a ~ 0.2 our theory yields a 
reasonable quantitative description of the effect. The ex- 
tension of our work for even stronger interaction and/or 
vanishing disorder will be published separatelyi^. 

3. Symmetric setup 

Another possibility^^, is to align the Dirac points in the 
two layers using a combination of gates and then apply a 
voltage V between the two layers, inducing same number 
of electrons in one layer and holes in another (such that 
ni = 77.2 = n and fii = — /Z2 = n)- In this case the results 
of Sec. He] and WJ\ apply. 

The inter-layer voltage V is related to the carrier den- 
sity by V = e^n/C + 2/i. The capacitance C is the only 
independently measurable coefficient in this relation and 
may be found by measuring the Hall coefficient at a large 
chemical potential. Then the V-dependence of any quan- 
tity can be directly translated into the density depen- 
dence (using Eq. (|65p to convert to n). 



V. CONCLUSIONS 

We have presented the perturbative theory of the 
Coulomb drag in ballistic graphene-based double-layer 
structures. Our theory is applicable to the wide range 
of temperatures and carrier densities, but is subject to 
the condition In addition we have limited our dis- 
cussion to the experimental^ condition ([3]), the former is 
necessary to justify the theoretical approach that we've 
adopted in this paper. As shown in Sec. lIVi Eq. ([5]) allows 
us to simplify our calculations by disregarding the effects 
of plasmon modes, Dirac spectrum renormalization, and 
energy dependent impurity scattering time. Eq. ^ also 
justifies our assumption that impurity scattering domi- 
nates the transport properties in the system. 

The main results of this paper can be summarized as 
follows. Qualitatively, poifJ^/T) has the same shape in 
all parameter regimes: in the vicinity of the Dirac point 
Pd /T'^, ^ ^ T the drag reaches its maximum and 
then decays at /i 3> T. This decay occurs over a wide 
region of /i where po cannot be described by a single 
power law. In particular, we have analyzed three regimes: 

(i) in the limit a ^ and d — > the drag coefficient 
is given by Eqs. ([23, dM]), and see also Table H] 

(ii) for a — >■ 0, but finite d the drag coefficient acquires 
logarithmic corrections, see Fig. [31 (iii) for intermediate 
interaction strength and ^ T we describe the crossover 
between the logarithmic and the Fermi-liquid behavior, 
see Eq. (|44ep . The latter occurs only at the largest values 
of /J,, such that xd ^ 1. Thus our theory describes all 
qualitative features observed in Ref. d. 

Formally our results are applicable in the limit of weak 
interaction. The actual value of a in physical graphene is 
still the subject of a debate. Recent experiments^!^ sug- 
gest that at experimentally relevant temperatures the ef- 
fective (or renormalized) interaction parameter is rather 
small. In addition, if one takes into account dielectric 
properties of the substrate and/or the insulating layer 
between the two graphene sheets in the double-layer 
device^^, then the effective value of a will be even smaller. 

For ultra-clean graphene, where transport is domi- 
nated by electron-electron interaction, our theory should 
be generalized for stronger interaction. In this case also 
the single-layer conductivity becomes non-trivial. In our 
opinion, the most adequate method for such calculations 
is the method of the kinetic equatio n^^'^^i^^ . Our work 
in this direction will be reported elsewherei^. 
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Appendix A: Non-linear susceptibility in graphene 



Here we derive the non-linear susceptibility (|14p and 
consider a few limiting cases. 

The general definition of the non-linear susceptibility is 
given by Ec^. (|lip . which we repeat here for convenience: 



I- Jdri j dr2T{uj-ri,r2)V{ri)V{r2). 



(Al) 



Furthermore, Eqs. p2|) and (|13l) express the non-linear 
susceptibility of a disordered conductor in terms of exact 
Green's functions of the system (for detailed derivation 
see Refs.ili3): 



de 
2^ 



( tanh ^J, tanh — — ) 7i2(e; ^) 



2T 



2T 



tanh 



2T 



e — cj — /i , 
tanh — ) (e; -co) 



, (A2) 



The factor iV = 4 reflects the spin and valley degeneracy, 
and 



k,k' 



1 



(A8) 



The current operator— should also be written in the ba- 
sis of the eigenstates where, unlike the original Bloch 
basis^^, the current operator depends on the direction of 
the quasi-particle momentum 



= 2eviynk, rik = k/k. 



(A9) 



The factor of 2 in Eq. (|A9[) appears due to the absence 
of backscattering in graphene: the transport time Ttr is 
twice the scattering time^^. 

Finally, in the ballistic regime 1/t and the Green's 
functions (IA6D can be written in the form 



G«(e;fc)G;^(e;fc) « 2nT{e)6{e - E,ik) 



(AlOa) 



7i2(e; ^) = [gU^- + c.) - GU^. + Gi (e) JgG^i (e) . 

(A3) 

In contrast to the usual Fermi Liquid calculation^ we 
shift the chemical potential from the Green's functions 
into the distribution functions. 

Averaging over disorder restores translational invari- 
ance. Moreover, in ballistic regime the Green's functions 
can be averaged independently. However, in graphene 
the eigenfunctions of the Dirac Hamiltonian are not plane 
waves. Focusing on a given valley and spin projection, 
we can write down the Dirac Hamiltonian as 



n^vk 



e*"^' 



cos ipk 



sm ipk 



(A4) 

The electron field operator in the basis of eigenstates can 
be written as 



(A5) 



where ^ = ± is the band index and the spinor [as well as 
the Hamiltonian (jA4l) ] is written in the sublattice space. 
In the basis of the eigenstates the (disorder-averaged) 
Green's functions are diagonal: 



G«(e,fc) = 



1 



e- E^{k)+i/2T{k)' 



Ey{k) = vvk. (A6) 



The impurity scattering time T(fe) may in general depend 
on momentum of the scattering states (see Appendix IB|) . 
Now we can write the triangular vertex 7 in the form: 



7(e;'^,Q) 



(A7) 



xImG« (e -I- c.; fc + g)G^ (e; k) 3Jg^ (e; k) . 



ImGf?(e; k) « -nSie - E^{k) 



(AlOb) 



Here we have replaced the momentum dependence of the 
scattering time by the energy dependence given the S- 
function approximation to the Green's functions. Note, 
that since Eq. (jA7p does not contain any energy inte- 
gration, this dependence plays no role in the triangular 
vertex 7, which in ballistic regime takes the form 

^{e;u;,q) = -NevT{e)J2'^ J rf^A: |a^';;^|'(A11) 
X(5(e — vvk)5 (e + cj - v'v\k + q\) . 



Using the (5-functions in Eq. (lAlip we notice, that in 
Eq. (|Alip the momenta satisfy 



(A12) 



k'^e'/v', {k + qy = {e + u;Y/v', 

and therefore 

_ (e + w)2 - - w^g^ _ + 2ecj - w^q^ , . 

Now we can replace the momentum dependence of the 
vertices A^'^_|_q by the frequency dependence: 

/ 2 



k,k-^-q 



= 1 



9 2 2 

uj — V q 



Then the triangular vertex (jAllI) becomes 



7(e;w,q) = —NevT{e) 



9 22 



(A14) 



(A15) 



g — ^ J (fk nkS{e — vvk)5 (e + a; — u'v\k + q\) . 

(A16) 
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Clearly, the direction of the vector g coincides with the f ^ sf , i \i , w 

J. J r 5o = > V I —;— o(e ~ vvk)b ie + Lo - V v\k + q\) . 

direction ot q: J k 

(gq) (A18) 

g = Aq =^ g = Q — , The remaining mtegration is straightforward and we find 

g2 



/kq 
(fk—6{e — vvk)5 (e + w — i''v\k + q\) . = —= 

Here we can again use Eq. (jA13|) and therefore 



4|e + w|0o(e,w,g) ^^^g^ 



g = q (A^l'^) where (6'(a;) is the Heaviside ^-function) 

I 



6'o(e,a;,9) = 6l(uq - |2e + a;|) [6l(-w - wg) - e[uj - vq)] + e{vq - \uj\) [e{2e + uj-vq) - e{-vq-2e-u)] . (A20) 
Now the triangular vertex gamma can be written using Eqs. (jAlSI) . (|A17|) . (|A19[) . and (|A20[) in the form ([T4| 



7 e,w,<7 =-q-Afe«r e „ ^ 2 \ ' 5^ sgn e + ^ 0^ e, ^, g . A21 



The function 0o(£j "^j l) is antisymmetric under the simultaneous change of sign of both frequencies 

%{>^,^,q) = -So{-e,~uj,q). 

Therefore the triangular vertex 7 as a whole is symmetric under the simultaneous change of sign of all variables: 

7(e,c^,q) =7(-e,-w,-q). (A22) 
Using this property in Eq. (|A2p we find the expression for the non-linear susceptibility in graphene: 

Nevq 



2tt 



g{uj,q;fj.), (A23a) 



g(w,<7;M) = J deTie)Iie,io)Fie,oj;q)eoie,oj,q), (A23b) 



/(e, ^) = tanh - tanh - tanh ^ + tanh . ( A23d) 



The expression then follows after a change of vari- 
ables indicated in the text preceding Eq. ([T6)) and explic- 
itly resolving the integration limits given by Eq. (|A20l) . 



Comparing Eq. (jA23p with the standard Fermi liquid 
result of Ref. [1, we should recall that within the Fermi 
liquid theory the Fermi energy is the largest energy scale 
in the problem. In order to compare Eq. (|A22I) to the 



corresponding results of Ref. [sj, we consider the limit 

w, ^ e ^ Ep- 

In this limit the leading contribution to the functions in 
the integrand in Eq. (jA23|) is given by 



F{uj,vq < e) 



we 
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FIG. 6: [Color online] A numerical comparison of the approx- 
imations (|29() and (jA25|l . The former is shown by the red 
solid line, while the latter - by the blue dashed line. The solid 
green line is the numerical evaluation of Eq. (IA23p . 



/2(cj«e)« — 



e{e) 



cosh 



2 £ 



££ 

2T 



cosh' 



2 e 



£i 

2T 



Thus all three terms F, I2, and in the integrand are 
odd functions of e and as a result the leading contribu- 
tion to non-linear susceptibility vanishesi^ii^. The sub- 
leading contribution (stemming from the second term in 
the function F) yields for vq ^ uj 



(A24) 



in agreement with Ref. 's'. Same result is given also by 
Eq. ([29]) in the limit /i > T for > w. 

The result (|A24l) was obtained assuming that typical 
values of momentum q are small compared to the Fermi 
momentum. This is justified for large inter-layer spacing 
>cd 1, which was assumed in Ref. S However, in 
graphene-based samples of Ref. |^ the inter-layer spacing 
is rather small and xd ~ 1. In this case we can no longer 
assume the momentum q to be small. Expanding the 
integrand in Eq. (|A23|) in the limit uj <^ vq,e ~ Ep, we 
find 



F{uj < vq,e) 



1 

2^ 



1 - 



1^' 



while the expansion of I2 remains the same as above. 
Now the frequency integration yields 



5(^«9)«-4-Jl-^. 



This expression approximates the non-linear susceptibil- 
ity in the region uj <^ vq < 2Ep. Comparing Eq. (|A25p 
to Eq. (|29|) we note that both approximations work well 
in that intermediate region, see Fig. [5] At the same time, 
Eq. ((29)) accounts better for the behavior at g ~ w and 
also allows for momenta larger than 2kF- Note, that 
the non-linear susceptibility (jA23[) is real, the fact that 
Eq. (jA25p yields imaginary values for q > 2kF is the 
artifact of its approximate derivation. 



Appendix B: Kinetic equation approach! to drag in 
graphene 

In this Appendix we derive the general expression for 
drag conductivity ao in the framework of the kinetic 
equation (this is justified by requiring large single-layer 
conductivity cri(2) ^ e^j which is valid for /i ^ T or 
/i » r). By solving the two coupled equations for the 
distribution functions of two graphene layers we will re- 
produce the result ^ with the nonlinear susceptibility 
T{uj,q) given by Eqs. ^ and pi)) . 

In Appendix |^ we have characterized the eigenstates 
()A5p of the massless Dirac Hamiltonian H ~ vcrk [see 
also Eq. ()A4p ] by the value of momentum k and the dis- 
crete variable v = ±1 indexing conduction and valence 
bands. In this representation, the electron energy and ve- 
locity are Ey{k) = vvk and v = vvek (where — k/k is 
the unit vector pointing in the direction of momentum) . 

Now for the purposes of deriving the kinetic equation, 
we find it more convenient to label the eigenstates by 
their energy e and the unit vector = v/v . The particle 
momentum is then k = evEiy{k)/v and the eigenstates 
are normalized as follows 



\e\dede-^ 

|e,e„)(e,e„| = 1. 



(27rw)2 



(Bl) 



In the lowest order of the perturbation theory we ne- 
glect electron-electron interaction within each layer and 
disregard the back action of the drag current in the pas- 
sive layer onto the distribution function in the active 
layer. As a result the kinetic equation for the active layer 
is effectively decoupled and has the form 



eEv 



dfa (fa) - fa 



de 



Taie) 



(B2) 



vq 



Here E is the applied electric field and (e) is the trans- 
port time due to disorder scattering. Index a refers to 
the active layer. The distribution function depends on e 
and e^ ; angular brackets denote averaging with respect 
to the velocity direction. Within linear response we sub- 
stitute the equilibrium distribution function in the 
(A25) left-hand side of Eq. ()B2I) and find the following solution 
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fa = - ra{e)^eEv = + r.(.)/(°) (l - ^ 



(0) 



eEv 



(B3) 



This is equation is written for the case of "dirty" graphene ^ '''eJ' ■ The opposite limit of clean graphene will be 
discussed in Ref. [H. 

Consider now the passive layer. We denote the corresponding distribution function /{, and include the collision term 
describing inter-layer scattering. The second kinetic equation has the form 



0: 



ifb) - h 



n 



+ J2 «^(«' a', b') k/^(l - /„)(! - h) - U,{1 - Q{l - fl) 



(B4) 



Here w{a^ b; a', b') is the probability of scattering (a', b') i— ^ (a, &), indices a, b, a' , and b' label incoming and scattered 
states in both layers. Summation over these states is carried out according to their normalization (|B1I) . 
The drag current can now be expressed as 



= eY,Vbfb^e J2 ^bVb w{a,b; a',b')[f'jUl - /a)(l - h) - fah{l - f'a){l - K) 



(B5) 



Now we substitute equilibrium distributions in layer 6, and the result (jB3p for /„ and Keeping only the terms linear 
in the external field E and using the momentum conservation law, we express the drag current as jg = a'^^EP. Using 
the time-reversal invariancc of the scattering probability, w{a,b; a',b') = w{a',b'; a,b), we represent the resulting 
drag conductivity in the symmetric form 



e 

2T 



(^fc^fc ^ nVb)a{T'aV'a " Ta^'a)/3 w{a, b] a', b') /^/^(l - - /fc). 



(B6) 



,6.b' 



Here all distribution functions are taken at thermal equilibrium and the superscripts are suppressed for brevity. Each 
of the four scattering times entering the above equation is taken at the corresponding energy. 
The transition probability w{a, 6; a', 6') can be written with the help of the Fermi golden rule: 



w{a,b- a',b') = \{aMU\a',b')\\2^)H{ea + eb^e[^-e'^,)5{K + kb-K-K), 



(B7) 



where the matrix element of the inter-layer interaction includes the Dirac factors 
velocities: 

I (a, 6|t/|a', &') r = \U{ka fcl) r , , 



re-written in terms of the 



(B8) 



With this matrix element, we can separate the quantities related to layers a and b in the expression for drag conduc- 
tivity (IB6|) . This allows us to represent it in the form of Eq. ([7]): 



ST J (27r)3 sinh^ 



(B9) 



2T 



Ta{q,u:) = {e-l^ - 1) ^ /^(l - /,)(t>:, - TaVa)^-^-^^ {2nf5{ea - 4 + ^) 5{K - K + q) (BIO) 



and the same formula for r^. Note the symmetry relations T{—q, —uj) = — r(<7, — w) = r(q, w). 

Let us now evaluate the expression (jBlOp . Using the energy- velocity basis and resolving the energy delta function, 
we represent T as an integral over e and over two velocity directions and . With equilibrium Fermi distribution 
functions, this yields 



T{q,Lo) 



Sirvq^ 



de\e{e + uj)\ 



tanh 



e + Lu ~ fi 
2T 



tanh 



e-/i 



2T 



J{e,e + uj,q), 



(Bll) 



J(e, e', l) ^ J de'^ {r'qe'^ - Tqe^){l + e^e'^)S{ee^ - e'e^ + vq) 



(B12) 
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The two-dimensional delta function in the latter integral fixes both e„ and e'^ . We substitute q from the argument of 
the delta function into the rest of the integrand and then average the delta function over directions of q. After such 
an averaging the integrand depends only on the angle between e„ and . 

J(e,e',g) = J [r'e' + re - (r'e + re') cos </>](!+ cos (^^J + e'^ - 2ee' cos(j)-vq^ = J(e,e',g) + J(e',e,q), 

(B13a) 

Ji^'^'l) - -(^) y .2g2_(,_,,)2 - (B13b) 

This result should be treated as zero if the argument of the square root is negative [this fact was previously expressed 
in terms of the additional factor 0o(e, w, q)]. Taking advantage of the symmetry of Eq. (jB13l) . we recast T in the form 
of Eqs. (in]) and ([M]). In particular, we identify the function (|B13|) with Eq. as 

7(e,^,9) = |e(e + ^^)l^(e,e + w,(?), 

where we have multiplied the result by N to account for the spin and valley degeneracy. Note, that in this Appendix 
T stands for the transport scattering time, unlike Eq. (jl4p . where r is just the mean free time. In graphene these two 
quantities differ by a factor of two2^. 



Appendix C: Polarization operator in graphene 



In the basis of exact eigenstates we can use the standard expression for the polarization operator, including the 
vertices Aj^'^, and summing over the two bands: 



[ J^W^y i2 nf(fc)-nf(fc + g) 

n {u;,q)--4^j :;-r^;(kTqH^Ak)Tl^ 



(Cl) 



[the prefactor of 4 is due to spin and valley degeneracy; the overall sign is chosen in such a way that the static 
polarization operator at 5 = yields the density of states (|63)) ]. Here np{k) stands for the Fermi distribution. 

The polarization operator was calculated in detail in Ref. 14. Nevertheless, we will add some details in order to make 
the paper self-contained. The complete expression for the polarization operator might also be useful for numerical 
computations. 

In order to simplify the expression for the polarization operator we multiply Eq. (|C1|) by the integral of a (5-function, 
which is unity: 



Now we can use the (5-functions to express the integrand in Eq. (|Cip in terms of e^. After that the momentum integral 
will only contain the two (5-functions and can be evaluated analytically similarly to how it was done in Appendix |X] 
for the non-linear susceptibility. Then the polarization operator takes the form 

n^ = - / nF eij-TiF £2 ■ -— F ei,e2), C3) 

J £162 UJ — €2 + €1 + 11] 



where 



^ /■ (Pk 

F{ei,e2) = Z2J " ^ - i^vlk + q\) . 



(C4) 



The calculation of this function can be performed along the lines of Appendix El The result is 



F{e„e2) = 'i'i ^2 2ir2 . ^ ^i^i^e^), (C5) 
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where 



e(ei, £2) = 0{ei > 0)9 [(ei - vq)^ < ej < (ei + vq)^] - 0{ei < 0)9 [(ei + vq)^ < < (ei - vq)^] 



(C6) 



The ^-functions are the result of imposing the condition that the cosine of the angle between the two momenta in 
Eq. (jA13[) is less than unity. In other words, the expression under the square root in Eq. (IC5[) has to be positive (and 
thus -F'(ei,e2) is a real function). 

The resulting expression can be simplified by the series of simple transformations: (i) change the sign of in the 
second term in Eq. (jC6|) : (ii) resolve the 0- functions in order to identify the integration limits; (iii) introduce the sum 
and difference 



Zi — Si + 62, Z-i — t\ — t2i 



(iv) introduce the dimensionless variables 



^ 2T' 2T T 



As a result we arrive at the following expression 
1 1 



n 



R 



dzidz2 



z,^{l-zl){l-zl) 



(^r^ - 1) 



Q 



Q 



Z2Q + W + ir] Z2Q -W - iri 

Q 



+ a-4) 



Jl(Zl ^,Z2,Xi) 

Q 



(C7) 



Zj^^Q + W + if] z^ 



W ~ irj 



J2{z^ ,Z2,Xi) 



T I \ i- ,i.Zi+Z2)Q + X ,{Z1+Z2)Q-X ,{Z1-Z2)Q + X (zi - Z2)Q - X 

Ji{2){zi, Z2,x) — tanh h tanh =F tanh =F tanh (C8) 

In particular, it is instructive to further simplify the imaginary part of the polarization operator: 



lraU''{iu,q) = [9{\W\ > Q)PiiW,Q) + 9{\W\ < Q)P2{W,Q)] ; 



Attv 



(C9a) 



1 

P^iW,Q) ^ -^^^= J dzVT^ hiz;Q,W,x) 



(C9b) 



P2{W,Q) = 9i\W\<Q) 



sgnW 



dzy/z^ - 1 /i(z;Q,W,x); 



(C9c) 



Ii{z]Q,W,x) — tanh ^ ^ 1_ ^^j^j^ Q + W — a; ^ ^^^^i £^ — ^ — — ~ tanh — ^ — W + x ^ (C9d) 



Comparing Eqs. (jC9p and ([T5)) we conclude, that despite clear similarity, these expression are not proportional 
to each other. Therefore the proportionality between the non-linear susceptibility and the imaginary part of the 
polarization operator mentioned in Ref. '3 is not a general theorem, but rather a property of the limiting cases 
considered in Ref. 

Having the full expression for the polarization operator it is straightforward to derive the well-known expressions: 



n^(;i = a; = 0; T < ug) w q/Av, 
n^(// = uj = 0;T:$>vq)M ATlii2/ (ttv^), 
U"{q < fc^; cj = T = 0) w 2kF/nv, 



(CIO) 
(Cll) 
(C12) 



Ren^((7<fcF;T = 0) 



2kF 

nv 



- vq 



Imn^(g< kF\T = 0) 



2kF Lo 9{\uj\ < vq) 

TTV y/v^q^ - Ul^ 



(C13) 
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Appendix D: Numerical evaluation of the drag coefficient 



In this Appendix we show the resuhs of the numerical evaluation of the drag coefficient using Eqs. ([TJ, and 
PU)) . The interaction propagator (jH]) was calculated using the polarization operator calculated in Appendix [Cl in the 
absence of disorder, see Eq. (jC7p . The non-linear susceptibility was evaluated using Eqs. (|16l) . The particle density 
was found from Eq. (j65p . The particular values of the inter- layer spacing d, temperature T, and the interaction 
parameter a were chosen to resemble possible realizations of the drag measurement in graphene-based devices^i^. 
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FIG. 7: Numerical estimate of the drag resistance. The left panel shows pb{^J-/T) for a = 0.3 (see Refl25h and various values 
of d (see Table lll| . The right panel shows pd(i-l/T) for d — 2Td/v and various values of a. 
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FIG. 8: The data from Fig. [7] in the linear scale. The chemical potential is converted to the carrier density, n, using Eq. (|65 



TABLE II: Values of the dimensionless parameter d = 2Td/hv for different inter-layer separations d and temperatures T. 



d = 4nm 



6 nm 



8 nm 



12 nm 



16 nm 



18 nm 



25 K 
50 K 
100 K 

200 K 



0.026 
0.052 
0.105 
0.209 



0.039 
0.078 
0.157 
0.314 



0.052 
0.105 
0.209 
0.418 



0.078 
0.157 
0.314 
0.628 



0.105 
0.209 
0.418 
0.837 



0.118 
0.235 
0.471 
0.942 
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The maximal values of are apparently reached for 
a ^ 0.5. The drag resistance at the peak is a non- 
monotonous function of a, since we are calculating the 
drag conductivity within the lowest-order perturbation 
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interaction propagator in order to describe screening ef- 
fects. The peak values of are achieved for carrier 
densities such that fi ^ T (only weakly depending on d). 
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